Ssid SH Nutrients
Ana Palacio
2020-04-22
General project set-up
# Get all libraries and sources required to run the script
library(plyr)
library(dplyr)
library(reshape2)
library(ggplot2)
library(ggthemes)
library(scales)
library(lme4)
library(lmerTest)
library(emmeans)
library(skimr)
# Default ggplot settings
Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))
ggthe_bw<-theme(plot.background=element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)+
theme_bw()1. DATA Exploration
1. Select qPCR info and define factors
# Import the data
qPCR.variables <- read.csv("Outputs/Ssid_SH_cell_ratio.csv", header = T)
# Remove blastate samples
qPCR.blastate<-qPCR.variables[((qPCR.variables$Date=="2018-02-05") |
(qPCR.variables$Date=="2018-03-08")), ]
qPCR.variables <- droplevels(qPCR.variables[!rownames(qPCR.variables) %in% rownames(qPCR.blastate), ])
# Variable types
str(qPCR.variables)## 'data.frame': 644 obs. of 17 variables:
## $ Treatment : Factor w/ 3 levels "A","N","N+P": 1 1 1 1 1 1 1 1 1 1 ...
## $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fragment : Factor w/ 162 levels "Ss_20-1","Ss_20-11",..: 27 30 31 48 49 52 55 58 97 100 ...
## $ Genotype : Factor w/ 7 levels "Ss_20","Ss_22",..: 2 2 2 2 3 3 3 3 5 5 ...
## $ Date : Factor w/ 5 levels "2017-11-15","2017-12-14",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Days : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Time_Point: Factor w/ 5 levels "T13","T19","T22",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Phase : Factor w/ 3 levels "Baseline","Heat",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Sample : Factor w/ 644 levels "Ss_20-1_T5","Ss_20-1_T9",..: 109 124 128 188 191 204 217 232 380 393 ...
## $ C.Ssid : num 0 0 0 0 0 ...
## $ D.Ssid : num 0.0266 0.034 0.0239 0.0281 1.0265 ...
## $ TotalSH : num 0.0266 0.034 0.0239 0.0281 1.0265 ...
## $ logC.SH : num NA NA NA NA NA ...
## $ logD.SH : num -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
## $ logSH : num -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
## $ D.Prp : num 1 1 1 1 1 ...
## $ Community : Factor w/ 3 levels "C1","C3","D": 3 3 3 3 3 3 3 3 3 3 ... qPCR.variables$Genotype<-factor(qPCR.variables$Genotype,
levels=c("Ss_22","Ss_23","Ss_27", "Ss_28",
"Ss_20", "Ss_24","Ss_30"))
qPCR.variables$DaysF<-as.factor(qPCR.variables$Days)
summary(qPCR.variables)## Treatment Replicate Fragment Genotype Date
## A :220 R1:333 Ss_20-17: 5 Ss_22:92 2017-11-15:140
## N :228 R2:311 Ss_20-19: 5 Ss_23:96 2017-12-14:154
## N+P:196 Ss_20-22: 5 Ss_27:96 2018-02-01:151
## Ss_20-26: 5 Ss_28:82 2018-02-23:103
## Ss_20-29: 5 Ss_20:97 2018-03-06: 96
## Ss_20-32: 5 Ss_24:85
## (Other) :614 Ss_30:96
## Days Time_Point Phase Sample
## Min. : 0.00 T13:151 Baseline :140 Ss_20-1_T5 : 1
## 1st Qu.: 29.00 T19:103 Heat :199 Ss_20-1_T9 : 1
## Median : 78.00 T22: 96 Nutrients:305 Ss_20-11_T13: 1
## Mean : 57.76 T5 :140 Ss_20-11_T19: 1
## 3rd Qu.:100.00 T9 :154 Ss_20-11_T22: 1
## Max. :111.00 Ss_20-11_T9 : 1
## (Other) :638
## C.Ssid D.Ssid TotalSH
## Min. :0.000000 Min. :0.000000 Min. :0.0000034
## 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.0093391
## Median :0.002373 Median :0.006613 Median :0.0313298
## Mean :0.046419 Mean :0.059660 Mean :0.1060796
## 3rd Qu.:0.027225 3rd Qu.:0.044028 3rd Qu.:0.1286860
## Max. :1.155796 Max. :1.394495 Max. :1.3944953
##
## logC.SH logD.SH logSH D.Prp
## Min. :-5.47264 Min. :-4.7761 Min. :-5.4726 Min. :0.0000
## 1st Qu.:-2.80684 1st Qu.:-2.2410 1st Qu.:-2.0297 1st Qu.:0.0000
## Median :-2.09749 Median :-1.6835 Median :-1.5040 Median :0.6792
## Mean :-2.08904 Mean :-1.7359 Mean :-1.5221 Mean :0.5242
## 3rd Qu.:-1.27217 3rd Qu.:-1.0631 3rd Qu.:-0.8905 3rd Qu.:1.0000
## Max. : 0.06288 Max. : 0.1444 Max. : 0.1444 Max. :1.0000
## NA's :184 NA's :205
## Community DaysF
## C1:132 0 :140
## C3:171 29 :154
## D :341 78 :151
## 100:103
## 111: 96
##
## skim(qPCR.variables)| Name | qPCR.variables |
| Number of rows | 644 |
| Number of columns | 18 |
| _______________________ | |
| Column type frequency: | |
| factor | 10 |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Treatment | 0 | 1 | FALSE | 3 | N: 228, A: 220, N+P: 196 |
| Replicate | 0 | 1 | FALSE | 2 | R1: 333, R2: 311 |
| Fragment | 0 | 1 | FALSE | 162 | Ss_: 5, Ss_: 5, Ss_: 5, Ss_: 5 |
| Genotype | 0 | 1 | FALSE | 7 | Ss_: 97, Ss_: 96, Ss_: 96, Ss_: 96 |
| Date | 0 | 1 | FALSE | 5 | 201: 154, 201: 151, 201: 140, 201: 103 |
| Time_Point | 0 | 1 | FALSE | 5 | T9: 154, T13: 151, T5: 140, T19: 103 |
| Phase | 0 | 1 | FALSE | 3 | Nut: 305, Hea: 199, Bas: 140 |
| Sample | 0 | 1 | FALSE | 644 | Ss_: 1, Ss_: 1, Ss_: 1, Ss_: 1 |
| Community | 0 | 1 | FALSE | 3 | D: 341, C3: 171, C1: 132 |
| DaysF | 0 | 1 | FALSE | 5 | 29: 154, 78: 151, 0: 140, 100: 103 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Days | 0 | 1.00 | 57.76 | 41.59 | 0.00 | 29.00 | 78.00 | 100.00 | 111.00 | ▆▆▁▆▇ |
| C.Ssid | 0 | 1.00 | 0.05 | 0.12 | 0.00 | 0.00 | 0.00 | 0.03 | 1.16 | ▇▁▁▁▁ |
| D.Ssid | 0 | 1.00 | 0.06 | 0.14 | 0.00 | 0.00 | 0.01 | 0.04 | 1.39 | ▇▁▁▁▁ |
| TotalSH | 0 | 1.00 | 0.11 | 0.17 | 0.00 | 0.01 | 0.03 | 0.13 | 1.39 | ▇▁▁▁▁ |
| logC.SH | 184 | 0.71 | -2.09 | 1.08 | -5.47 | -2.81 | -2.10 | -1.27 | 0.06 | ▁▃▇▇▅ |
| logD.SH | 205 | 0.68 | -1.74 | 0.94 | -4.78 | -2.24 | -1.68 | -1.06 | 0.14 | ▁▂▆▇▃ |
| logSH | 0 | 1.00 | -1.52 | 0.83 | -5.47 | -2.03 | -1.50 | -0.89 | 0.14 | ▁▁▃▇▅ |
| D.Prp | 0 | 1.00 | 0.52 | 0.47 | 0.00 | 0.00 | 0.68 | 1.00 | 1.00 | ▇▁▁▁▇ |
2. Exploratory graphs
- Histograms
HistoL_SH<-qplot(logSH, data=qPCR.variables, binwidth=0.15)
HistoL_SH + facet_grid(Treatment~Date)
ggplot(qPCR.variables, aes(logSH, fill = Treatment , colour = Treatment)) +
geom_density(alpha = 0.1) + facet_wrap(~Date) + ggthe_bw
- Log SH
logSHGenotype <- ggplot(qPCR.variables, aes(Days, logSH)) +
geom_line(aes(colour=Fragment))+geom_point(aes(shape=factor(Community), colour=factor(Community)))+
# geom_jitter(aes(colour=factor(Replicate))) +
# stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
#stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
# stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) +
facet_grid(Treatment~Genotype) +
ggthe_bw +theme(legend.position = "none" )
logSHGenotype + ylab("Relative log10 (S:H)") + xlab("Treatment") +
theme(axis.title.y=element_text(size=12), legend.position="none")
logSHTreatment <- ggplot(qPCR.variables, aes(Date, logSH, colour=factor(Treatment))) +
# geom_jitter(aes(colour=factor(Treatment))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
stat_summary(fun.y=mean, geom="line", size =1) +
theme_bw()
logSHTreatment + facet_wrap(~Genotype)
logSH_Replicate<- ggplot(qPCR.variables, aes (Days, logSH,
colour=factor(Replicate))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
stat_summary(fun.y=mean, geom="line") + facet_grid (~Treatment) + ggthe_bw
logSH_Replicate 
logSH_Replicate+ facet_grid(~Community)
logSHTreatment <- ggplot (qPCR.variables, aes(Days, logSH, colour=Treatment)) +
ggtitle("A.") + # geom_point(alpha=0.3) +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
#stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
# position = position_dodge(width = 1)) +
stat_summary(fun.y=mean, geom="line", position = position_dodge(width = 5)) +
scale_y_continuous(breaks = seq(-5, 0.5, 0.5),
expand = c(0,0),
name=("log 10 (Total S/H cell ratio")) +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 15),
expand = c(0, 0)) +
annotate("segment", x = 2, xend = 91, y = -3.5, yend = -3.5,
colour = "gray90")+
annotate("segment", x = 79, xend = 90, y = -3.5, yend = -3,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y = -3, yend = -3,
colour = "gray90")
logSHTreatment
#logSHTreatment + facet_wrap(~Genotype)
logSHTreatment + facet_wrap(~Community)
- Log C
logC_Treatment <- ggplot(qPCR.variables,
aes(Days, logC.SH, colour=factor(Fragment))) +
geom_line()+
#geom_jitter(aes(colour=factor(Replicate))) +
#stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
#stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) +
ggthe_bw + theme(legend.position = "none" )
logC_Treatment
logC_Treatment + facet_grid(Genotype~Treatment)
logC_Treatment + facet_grid(Community~Treatment)
- Log D
logD_genotype <- ggplot(qPCR.variables,
aes(Days, logD.SH, colour=factor(Fragment))) +
geom_line()+
#geom_jitter(aes(colour=factor(Replicate))) +
#stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
#stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) +
ggthe_bw + theme(legend.position = "none" )
logD_genotype
logD_genotype + facet_grid(Genotype~Treatment)
logD_Treat <- ggplot(qPCR.variables,
aes(Days, logD.SH, colour=factor(Treatment))) +
#geom_jitter(aes(colour=factor(Replicate))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) +
ggthe_bw + theme(legend.position = "none" )
logD_Treat
logD_Treat + facet_grid(Community~Treatment)
*D proportion
qPCR.variables$D.Prp[qPCR.variables$Sample=="Ss_22-18_T22"]<-NA
qPCR.variables$Community[qPCR.variables$Sample=="Ss_22-18_T22"]<-"D"
qPCR.variables$D.Prp[qPCR.variables$Sample=="Ss_23-80_T19"]<-NA
qPCR.variables$Community[qPCR.variables$Sample=="Ss_23-80_T19"]<-"D"
D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
ggthe_bw + theme(legend.position="none") +
geom_jitter(aes(colour=factor(Fragment)), alpha=0.5, size=1) +
geom_line() + facet_grid (Genotype~Treatment) +
annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
colour = "gray90")+
annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
colour = "gray90")+
annotate("text", x = 8, y = 0.08, label = "(BL)", size=3, colour="gray")+
annotate("text", x = 45, y = 0.08, label = "(C)", size=3, colour="gray")+
annotate("text", x = 100, y = 0.08, label = "(H)", size=3, colour="gray")+
scale_y_continuous(breaks = seq(0, 1, 0.3),
name=("Durusdinium proportion (D/H)/(S/H)")) +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 45),
expand = c(0, 0))
D.PTreatment
ggsave(file="Outputs/D.P_Treatment_Genotype.svg", plot=D.PTreatment, width=4.5, height=8)2. Subset data
# Removing weird points and colonies not used
qPCR.variables_2<-subset(qPCR.variables, Date!="2017-12-14")
qPCR.variables_2<-subset(qPCR.variables_2, TotalSH<0.8)
SH.0<-subset(qPCR.variables_2, Days<2)
SH.C<-subset(qPCR.variables_2, Days<79)
SH.C<-subset(SH.C, Days>77)
SH.H<-subset(qPCR.variables_2, Days>80)
logSH_Genotype<- ggplot(qPCR.variables_2, aes (Days, logSH,
colour=factor(Fragment))) +
theme_bw() + geom_line()+ facet_grid (Genotype~Treatment)#+
logSH_Genotype + theme(legend.position = "none" ) + geom_point()
Log S/H graph
logSHTreatment <- ggplot (qPCR.variables_2, aes(Days, logSH, colour=Treatment)) + theme_bw() + ggtitle("A.")+ Fill.colour+ ggthe_bw+
theme(plot.background=element_blank(),
panel.border = element_blank(),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position="bottom",
strip.background = element_rect(fill="white")) +
#geom_jitter(aes(colour=factor(Replicate))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
#stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
# position = position_dodge(width = 5)) +
stat_summary(fun.y=mean, geom="line") +
annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
colour = "gray90", linetype=1)+
annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
colour = "gray90", linetype=1)+
annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
colour = "gray90", linetype=1)+
annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
scale_y_continuous(breaks = seq(-5, 0.3, 0.5),
expand = c(0,0),
name=("Symbiodiniaceae / S. siderea (S/H) cell ratio")) +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 15),
expand = c(0, 0))
logSHTreatment 
logSHTreatment + facet_grid(~Community)
#ggsave(file="5.3A_AH.svg", plot=logSHTreatment, width=3.0, height=6)
#ggsave(file="5.3A_AH_NoGenotype.svg", plot=logSHTreatment, width=3.0, height=3)Log C/H graph
CHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logC.SH, colour=Treatment)) + theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
colour = "gray90", linetype=1)+
annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
colour = "gray90", linetype=1)+
annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
colour = "gray90", linetype=1)+
annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
#geom_jitter(aes(colour=factor(Replicate))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
#stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
# position = position_dodge(width = 5)) +
stat_summary(fun.y=mean, geom="line") +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 15),
expand = c(0, 0))
CHTreatment_2 
#CHTreatment_2 + facet_grid (Genotype~Community)
CHTreatment_2 + facet_grid(~Community)
#ggsave(file="5.3A_AH_2.svg", plot=SHTreatment_2, width=3.0, height=6)
#ggsave(file="5.3C_SH_NoGenotype_2.svg", plot=SHTreatment_2, width=3.0, height=3)Log D/H graph
DHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logD.SH, colour=Treatment)) + theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
colour = "gray90", linetype=1)+
annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
colour = "gray90", linetype=1)+
annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
colour = "gray90", linetype=1)+
annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
#geom_jitter(aes(colour=factor(Replicate))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
#stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
# position = position_dodge(width = 5)) +
stat_summary(fun.y=mean, geom="line") +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 15),
expand = c(0, 0))
DHTreatment_2 + facet_grid (Genotype~Community)
DHTreatment_2 + facet_grid(~Community)
D Proportion
D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
ggthe_bw + theme(legend.position="none") +
geom_jitter(aes(shape=factor(Treatment)), alpha=0.5) +
geom_line()+
annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
colour = "gray90")+
#annotate("text", x = 45, y = 0.2, label = "Nutrients", size=3)+
annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
colour = "gray90")+
#annotate("text", x = 99, y = 0.3, label = "Heat", size=3)+
scale_y_continuous(breaks = seq(0, 1, 0.3),
name=("Durusdinium proportion (D/H)/(S/H)")) +
scale_x_continuous(name="Days in the experiment",
limits = c(-2,113),
breaks = seq(0, 110, 45),
expand = c(0, 0))
D.PTreatment + facet_grid (Genotype~Replicate)
D.PTreatment
3. GLM models
Baseline
LME_BaseLine<-lmer(logSH ~ Treatment + Community + (1|Genotype) + (1|Replicate),
data=SH.0)
step(LME_BaseLine)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -75.016 166.03
## (1 | Genotype) 0 7 -82.009 178.02 13.986 1 0.0001842 ***
## (1 | Replicate) 0 7 -109.286 232.57 68.539 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 1 0.59627 0.29813 2 126.566 2.0641 0.1312
## Community 2 0.79889 0.39945 2 6.535 2.7209 0.1381
##
## Model found:
## logSH ~ (1 | Genotype) + (1 | Replicate) anova(LME_BaseLine) # Treatments is not significant ranova(LME_BaseLine) # Treatments is significant! LME_BaseLine<-lmer(logSH ~ Genotype + (1|Replicate),
data=SH.0)
step(LME_BaseLine)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 -72.767 163.53
## (1 | Replicate) 0 8 -106.895 229.79 68.256 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Genotype 0 8.6374 1.4396 6 129 9.7555 7.505e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Genotype + (1 | Replicate) anova(LME_BaseLine) # Treatments is not significant ranova(LME_BaseLine) # Treatments is significant!# Multicomp
Ssid.SH.C.emm<-emmeans(LME_BaseLine, ~Genotype)
#contrast(Ofav.SH.C.emm, "tukey")
Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
Ssid.SH.C_groups# Effect plot
plot(emmeans(LME_BaseLine, ~Genotype), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()
logSHTreatmentBL <- ggplot (SH.0, aes(Genotype, logSH, colour=Genotype)) +
#ggtitle("A.") +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position = position_dodge(width = 1)) +
scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio"))
logSHTreatmentBL
logSHTreatmentBL + facet_grid(~Community)
C: Effect of nutrient treatments at control temperature
LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate) + (1|Genotype),
data=SH.C)
step (LME_Ssid.C) # Remove Genotype## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -73.147 162.29
## (1 | Genotype) 1 7 -73.264 160.53 0.235 1 0.6278
## (1 | Replicate) 0 6 -105.884 223.77 65.240 1 6.63e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 1 0.7873 0.3936 2 144 2.8627 0.06037 .
## Community 0 9.4779 4.7390 2 146 33.6053 9.879e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Community + (1 | Replicate) anova (LME_Ssid.C) ranova (LME_Ssid.C)LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate),
data=SH.C)
step (LME_Ssid.C) # Remove replicate## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 -73.264 160.53
## (1 | Replicate) 0 6 -105.884 223.77 65.24 1 6.63e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 1 0.7873 0.3936 2 144 2.8627 0.06037 .
## Community 0 9.4779 4.7390 2 146 33.6053 9.879e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Community + (1 | Replicate) anova (LME_Ssid.C) ranova (LME_Ssid.C)# Multicomp
Ssid.SH.C.emm<-emmeans(LME_Ssid.C, ~Community)
#contrast(Ssid.SH.C.emm, "tukey")
Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
Ssid.SH.C_groups# Effect plot
plot(emmeans(LME_Ssid.C, ~Treatment), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()
logSHTreatmentC <- ggplot (SH.C, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) +
ggtitle("A.")+ Fill.colour +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position = position_dodge(width = 1)) +
scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio"))
logSHTreatmentC
logSHGenotypeC <- ggplot (SH.C, aes(Genotype, logSH, colour=Genotype)) +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position = position_dodge(width = 1)) +
scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (C/H) cell ratio"))
logSHGenotypeC
H: Effect of pre-exposure to nutrient treatments during heat challenge
LME_Ssid.H<-lmer(logSH ~ Treatment * Community * DaysF+ (1|Replicate) + (1|Genotype),
data=SH.H)
step (LME_Ssid.H) # Remove Replicate## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 21 -120.54 283.08
## (1 | Replicate) 1 20 -120.54 281.08 0.000 1 1
## (1 | Genotype) 0 19 -157.51 353.02 73.949 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:Community:DaysF 0 1.8668 0.4667 4 174.72 3.0217
## Pr(>F)
## Treatment:Community:DaysF 0.01928 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + Community + DaysF + (1 | Genotype) + Treatment:Community +
## Treatment:DaysF + Community:DaysF + Treatment:Community:DaysF anova (LME_Ssid.H) ranova (LME_Ssid.H)LME_Ssid.H<-lmer(logSH ~ Treatment * Community * DaysF+ (1|Genotype),
data=SH.H)
step (LME_Ssid.H) # Remove replicate## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 20 -120.54 281.08
## (1 | Genotype) 0 19 -157.51 353.02 73.949 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:Community:DaysF 0 1.8668 0.4667 4 174.72 3.0217
## Pr(>F)
## Treatment:Community:DaysF 0.01928 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment * Community * DaysF + (1 | Genotype) anova (LME_Ssid.H) ranova (LME_Ssid.H)# Multicomp
Ssid.SH.H.emm<-emmeans(LME_Ssid.H, ~Treatment*Community*DaysF)
#contrast(Ssid.SH.H.emm, "tukey")
Ssid.SH.H_groups<-cld(Ssid.SH.H.emm, by=NULL) # compact-letter display
Ssid.SH.H_groups<-Ssid.SH.H_groups[order(Ssid.SH.H_groups$Day,
Ssid.SH.H_groups$Community,
Ssid.SH.H_groups$Treatment), ]
Ssid.SH.H_groups# Effect plot
plot(emmeans(LME_Ssid.H, ~DaysF*Community|Treatment), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()
logSHTreatmentH <- ggplot (SH.H, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) +
ggtitle("A.")+ Fill.colour +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position = position_dodge(width = 1)) +
scale_y_continuous(breaks = seq(-3, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio"))
logSHTreatmentH
logSHTreatmentH + facet_grid(~Community)
logSHTreatmentH + facet_grid(DaysF~Community) #+ geom_jitter()
logSHTreatmentH <- ggplot (SH.H, aes(DaysF, logSH)) +
ggtitle("A.")+ Fill.colour +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position = position_dodge(width = 5), width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position = position_dodge(width = 1)) +
scale_y_continuous(breaks = seq(-3, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio"))
logSHTreatmentH
logSHTreatmentH + facet_grid(~Community)
logSHTreatmentH + facet_grid(Treatment~Community) #+ geom_jitter()
logSHGenotypeH <- ggplot (SH.H, aes(Community, logSH, colour=Genotype)) +
ggthe_bw + theme(legend.position="bottom") +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio"))
logSHGenotypeH
logSHGenotypeH + facet_grid(~DaysF)
All phases
- Dasys as a factor
LM_Ssid.qpCR<-lmer(logSH ~ Treatment * DaysF * Community + (1|Replicate) + (1|Genotype/Fragment),
data= qPCR.variables_2)
step(LM_Ssid.qpCR)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df
## <none> 40 -347.95 775.90
## (1 | Replicate) 1 39 -347.95 773.90 0.000 1
## (1 | Fragment:Genotype) 2 38 -347.95 771.90 0.000 1
## (1 | Genotype) 0 37 -377.52 829.03 59.127 1
## Pr(>Chisq)
## <none>
## (1 | Replicate) 1
## (1 | Fragment:Genotype) 1
## (1 | Genotype) 1.478e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:DaysF:Community 1 2.0375 0.1698 12 443.32 0.7963
## Treatment:Community 2 0.8016 0.2004 4 455.42 0.9449
## Treatment:DaysF 0 5.1227 0.8538 6 459.60 4.0277
## DaysF:Community 0 24.7264 4.1211 6 459.71 19.4411
## Pr(>F)
## Treatment:DaysF:Community 0.6545075
## Treatment:Community 0.4376783
## Treatment:DaysF 0.0006045 ***
## DaysF:Community < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF +
## DaysF:Community anova(LM_Ssid.qpCR) ranova(LM_Ssid.qpCR)LM_Ssid.qpCR<-lmer(logSH ~ Treatment + DaysF + Community +
(1 | Genotype) +
Treatment:DaysF + DaysF:Community,
data= qPCR.variables_2)
step(LM_Ssid.qpCR)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 -346.54 737.08
## (1 | Genotype) 0 21 -376.97 795.94 60.855 1 6.144e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment:DaysF 0 5.1227 0.8538 6 459.60 4.0277 0.0006045
## DaysF:Community 0 24.7264 4.1211 6 459.71 19.4411 < 2.2e-16
##
## Treatment:DaysF ***
## DaysF:Community ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF +
## DaysF:Community anova(LM_Ssid.qpCR) ranova(LM_Ssid.qpCR) summary(LM_Ssid.qpCR)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF +
## DaysF:Community
## Data: qPCR.variables_2
##
## REML criterion at convergence: 693.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3953 -0.6798 -0.0318 0.7153 3.3216
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.4236 0.6508
## Residual 0.2120 0.4604
## Number of obs: 486, groups: Genotype, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.89344 0.28407 8.84781 -6.665 9.97e-05 ***
## TreatmentN -0.07281 0.09363 459.23851 -0.778 0.437183
## TreatmentN+P -0.23487 0.09896 459.27507 -2.373 0.018035 *
## DaysF78 0.18335 0.14590 459.31051 1.257 0.209516
## DaysF100 -0.89065 0.15375 459.46923 -5.793 1.28e-08 ***
## DaysF111 -1.12847 0.15765 459.46253 -7.158 3.27e-12 ***
## CommunityC3 1.42388 0.22411 387.61081 6.353 5.91e-10 ***
## CommunityD 0.14326 0.14828 465.63001 0.966 0.334479
## TreatmentN:DaysF78 0.09799 0.13142 459.23892 0.746 0.456285
## TreatmentN+P:DaysF78 0.08542 0.13575 459.26579 0.629 0.529517
## TreatmentN:DaysF100 0.48243 0.14265 459.33134 3.382 0.000781 ***
## TreatmentN+P:DaysF100 0.44027 0.15202 459.29481 2.896 0.003957 **
## TreatmentN:DaysF111 0.09349 0.14708 460.66389 0.636 0.525338
## TreatmentN+P:DaysF111 0.46388 0.15488 459.43015 2.995 0.002891 **
## DaysF78:CommunityC3 0.35719 0.16195 459.31568 2.206 0.027907 *
## DaysF100:CommunityC3 -0.11623 0.17443 459.35517 -0.666 0.505507
## DaysF111:CommunityC3 -1.21190 0.19904 459.31114 -6.089 2.40e-09 ***
## DaysF78:CommunityD 0.10888 0.14638 459.32666 0.744 0.457370
## DaysF100:CommunityD 0.03883 0.15790 459.35328 0.246 0.805831
## DaysF111:CommunityD 0.26168 0.15986 460.69394 1.637 0.102321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Emmeans
Ssid.YII.emm<-emmeans(LM_Ssid.qpCR, ~Treatment*DaysF*Community)
#contrast(Ssid.YII.emm, "tukey")
# Pairwise comparisons
#pairs(Ssid.YII.emm) # same than contrast(Ssid.YII.emm, "tukey")
plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF), comparisons = TRUE) +
theme_bw() # Tukey comparission (do not trust CI to compare EMMs)
plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF|Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(Community~DaysF) +
theme_bw()
#CLD
SH_Multicomp<-cld(Ssid.YII.emm, by=NULL) # compact-letter display
SH_Multicomp<-SH_Multicomp[order(SH_Multicomp$DaysF,
SH_Multicomp$Community,
SH_Multicomp$Treatment), ]
#write.csv(SH_Multicomp, "Outputs/Ssid_SH_Multicomp.csv") # 2. Predict values:
pred_Ssid1 <- predict(LM_Ssid.qpCR, re.form = NA)
#3. Bootstrap CI:
Ss.boot1 <- bootMer(LM_Ssid.qpCR, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.boot1$t, 2, sd)
CI.lo_1 <- pred_Ssid1 - std.err*1.96
CI.hi_1 <- pred_Ssid1 + std.err*1.96
#Plot
Model_Ss_1b_plot<- ggplot(
qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
geom_line(aes(y = pred_Ssid1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
#scale_y_continuous(name=expression(~italic("Fv / Fm")),
# limits = c(0.5, 0.61),
# breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
#scale_x_continuous("Days in the experiment", limits = c(0, 78),
# breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1)
# stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw## List of 69
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.grid.major.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.grid.major.y : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.grid.minor.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.grid.minor.y : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 2.75pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 2.75pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 2.75pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 2.2pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 2.2pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 2.2pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.length : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.y : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing : 'unit' num 11pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'unit' num 1.2lines
## ..- attr(*, "valid.unit")= int 3
## ..- attr(*, "unit")= chr "lines"
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "valid.unit")= int 1
## ..- attr(*, "unit")= chr "cm"
## $ legend.box.spacing : 'unit' num 11pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'unit' num 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ size : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.ontop : logi FALSE
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4pt 4.4pt 4.4pt 4.4pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.switch.pad.grid : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.switch.pad.wrap : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUE Model_Ss_1b_plot + facet_wrap(~Community)
- Days as continous
LM_Ssid.qpCR2<-lmer(logSH ~ Treatment * Days * Community + (1|Replicate) + (1|Genotype/Fragment),
data= qPCR.variables_2)
step(LM_Ssid.qpCR2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 -568.24 1180.5
## (1 | Replicate) 1 21 -568.24 1178.5 0.00 1 1
## (1 | Fragment:Genotype) 2 20 -568.24 1176.5 0.00 1 1
## (1 | Genotype) 0 19 -589.20 1216.4 41.93 1 9.461e-11
##
## <none>
## (1 | Replicate)
## (1 | Fragment:Genotype)
## (1 | Genotype) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:Days:Community 1 1.7016 0.4254 4 461.15 0.9054
## Days:Community 2 1.7500 0.8750 2 465.80 1.8642
## Treatment:Community 3 3.2309 0.8077 4 467.27 1.7133
## Community 0 22.9355 11.4678 2 364.82 24.1860
## Treatment:Days 0 3.0544 1.5272 2 471.54 3.2209
## Pr(>F)
## Treatment:Days:Community 0.4606
## Days:Community 0.1562
## Treatment:Community 0.1458
## Community 1.369e-10 ***
## Treatment:Days 0.0408 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days anova(LM_Ssid.qpCR2) ranova(LM_Ssid.qpCR2)LM_Ssid.qpCR2<-lmer(logSH ~Treatment + Days + Community + (1|Genotype) + Treatment:Days,
data= qPCR.variables_2)
step(LM_Ssid.qpCR2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 10 -542.81 1105.6
## (1 | Genotype) 0 9 -564.04 1146.1 42.471 1 7.173e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Community 0 22.9355 11.4678 2 364.82 24.1860 1.369e-10
## Treatment:Days 0 3.0544 1.5272 2 471.54 3.2209 0.0408
##
## Community ***
## Treatment:Days *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days anova(LM_Ssid.qpCR2) ranova(LM_Ssid.qpCR2) summary(LM_Ssid.qpCR2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
## Data: qPCR.variables_2
##
## REML criterion at convergence: 1085.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5286 -0.6117 -0.1087 0.6507 2.3774
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.6803 0.8248
## Residual 0.4741 0.6886
## Number of obs: 486, groups: Genotype, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.853179 0.360321 8.344461 -5.143 0.000773 ***
## TreatmentN -0.096926 0.136250 471.094970 -0.711 0.477201
## TreatmentN+P -0.266358 0.143608 471.050405 -1.855 0.064257 .
## Days -0.008168 0.001209 471.227691 -6.757 4.18e-11 ***
## CommunityC3 1.818818 0.283439 290.421477 6.417 5.64e-10 ***
## CommunityD 0.262047 0.169203 453.518949 1.549 0.122148
## TreatmentN:Days 0.003229 0.001701 471.924639 1.898 0.058301 .
## TreatmentN+P:Days 0.004293 0.001799 471.108284 2.386 0.017420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnN TrtN+P Days CmmnC3 CmmntD TrtN:D
## TreatmentN -0.183
## TreatmntN+P -0.185 0.480
## Days -0.248 0.591 0.564
## CommunityC3 -0.371 -0.026 0.005 0.053
## CommunityD -0.382 -0.009 0.011 0.047 0.581
## TrtmntN:Dys 0.149 -0.834 -0.400 -0.706 0.063 -0.018
## TrtmntN+P:D 0.141 -0.399 -0.839 -0.669 0.021 0.020 0.477 # 2. Predict values:
pred_Ssid1 <- predict(LM_Ssid.qpCR2,re.form = NA)
#3. Bootstrap CI:
Ss.boot1 <- bootMer(LM_Ssid.qpCR2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.boot1$t, 2, sd)
CI.lo_1 <- pred_Ssid1 - std.err*1.96
CI.hi_1 <- pred_Ssid1 + std.err*1.96
#Plot
Model_Ss_1b_plot<- ggplot(
qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
geom_line(aes(y = pred_Ssid1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
#scale_y_continuous(name=expression(~italic("Fv / Fm")),
# limits = c(0.5, 0.61),
# breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
#scale_x_continuous("Days in the experiment", limits = c(0, 78),
# breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1) +
# stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_1b_plot + facet_wrap(Genotype~Community)
All phases, except baseline
- Dasys as a factor
qPCR.variables_3<-qPCR.variables_2[(qPCR.variables_2$DaysF!="0"), ]
Sy.Summary <- plyr::ddply (qPCR.variables_3, . (DaysF, Community, Treatment),
summarise,
Sy_mean = mean (TotalSH, na.rm = T),
Sy_sd = sd (TotalSH, na.rm = T))
Sy.Summarywrite.csv(Sy.Summary, "Outputs/meanSH.csv")Treatment<- ggplot(
qPCR.variables_3, aes(x = Days, y = logSH,
colour = Treatment, shape = Treatment)) +
Fill.colour+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(3) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(3),
linetype=1, alpha=1)+
stat_summary(fun.y=mean, geom="point", size =1,
position=position_dodge(3), alpha=0.8) +
ggthe_bw + facet_grid(~Community) +
scale_y_continuous(breaks = seq(-5.5, 0.5, 0.5),
expand = c(0.1,0.1),
name=("log 10 (Total S/H cell ratio)")) +
scale_x_continuous(name="Days in the experiment",
limits = c(75,113),
breaks = seq(75, 110, 15),
expand = c(0, 0)) +
theme(legend.position = c(.80, 0.24))+
annotate("text", x = 78, y = -5, label = "(C)", size=3, colour="gray")+
annotate("text", x = 100, y = -5, label = "(H)", size=3, colour="gray")
Treatment
#ggsave(file="Outputs/TotSH_Community.svg", plot=Treatment, width=4.5, height=3.5)LM_Ssid.qpCR<-lmer(logSH ~ Treatment * DaysF * Community + (1|Replicate) + (1|Genotype/Fragment),
data= qPCR.variables_3)
step(LM_Ssid.qpCR)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df
## <none> 31 -220.31 502.61
## (1 | Fragment:Genotype) 1 30 -220.31 500.61 0.000 1
## (1 | Replicate) 0 29 -237.04 532.09 33.477 1
## (1 | Genotype) 0 29 -249.84 557.68 59.068 1
## Pr(>Chisq)
## <none>
## (1 | Fragment:Genotype) 0.997
## (1 | Replicate) 7.210e-09 ***
## (1 | Genotype) 1.523e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:DaysF:Community 1 2.3236 0.2905 8 314.52 1.6881
## Treatment:Community 2 0.8550 0.2138 4 322.58 1.2212
## Treatment:DaysF 0 3.2886 0.8221 4 326.73 4.6834
## DaysF:Community 0 23.9600 5.9900 4 326.81 34.1226
## Pr(>F)
## Treatment:DaysF:Community 0.100414
## Treatment:Community 0.301634
## Treatment:DaysF 0.001092 **
## DaysF:Community < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 |
## Genotype) + Treatment:DaysF + DaysF:Community anova(LM_Ssid.qpCR) ranova(LM_Ssid.qpCR)LM_Ssid.qpCR<-lmer(logSH ~ Treatment + DaysF + Community +
(1 | Replicate) + (1 | Genotype) +
Treatment:DaysF + DaysF:Community, data= qPCR.variables_3)
step(LM_Ssid.qpCR)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 18 -222.50 480.99
## (1 | Replicate) 0 17 -237.29 508.58 29.595 1 5.325e-08 ***
## (1 | Genotype) 0 17 -251.70 537.40 58.408 1 2.130e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment:DaysF 0 3.2886 0.8221 4 326.73 4.6834 0.001092
## DaysF:Community 0 23.9600 5.9900 4 326.81 34.1226 < 2.2e-16
##
## Treatment:DaysF **
## DaysF:Community ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 |
## Genotype) + Treatment:DaysF + DaysF:Community anova(LM_Ssid.qpCR) ranova(LM_Ssid.qpCR) summary(LM_Ssid.qpCR)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 |
## Genotype) + Treatment:DaysF + DaysF:Community
## Data: qPCR.variables_3
##
## REML criterion at convergence: 445
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4207 -0.6488 0.1561 0.6873 2.2624
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.3836 0.6193
## Replicate (Intercept) 0.0364 0.1908
## Residual 0.1755 0.4190
## Number of obs: 349, groups: Genotype, 7; Replicate, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.646837 0.300750 8.678241 -5.476 0.000445
## TreatmentN 0.026826 0.083953 326.392615 0.320 0.749526
## TreatmentN+P -0.155045 0.084638 326.495151 -1.832 0.067885
## DaysF100 -1.047853 0.137310 326.459574 -7.631 2.58e-13
## DaysF111 -1.308796 0.141498 326.622442 -9.250 < 2e-16
## CommunityC3 1.660629 0.206771 291.644381 8.031 2.40e-14
## CommunityD 0.199342 0.138081 332.129562 1.444 0.149778
## TreatmentN:DaysF100 0.394931 0.129108 326.496928 3.059 0.002405
## TreatmentN+P:DaysF100 0.313005 0.134932 326.499993 2.320 0.020972
## TreatmentN:DaysF111 0.003397 0.133201 327.432319 0.026 0.979667
## TreatmentN+P:DaysF111 0.332339 0.137490 326.610740 2.417 0.016188
## DaysF100:CommunityC3 -0.492060 0.152909 326.388468 -3.218 0.001420
## DaysF111:CommunityC3 -1.576303 0.176948 326.557888 -8.908 < 2e-16
## DaysF100:CommunityD -0.126905 0.138844 326.531911 -0.914 0.361387
## DaysF111:CommunityD 0.140488 0.141114 327.593452 0.996 0.320197
##
## (Intercept) ***
## TreatmentN
## TreatmentN+P .
## DaysF100 ***
## DaysF111 ***
## CommunityC3 ***
## CommunityD
## TreatmentN:DaysF100 **
## TreatmentN+P:DaysF100 *
## TreatmentN:DaysF111
## TreatmentN+P:DaysF111 *
## DaysF100:CommunityC3 **
## DaysF111:CommunityC3 ***
## DaysF100:CommunityD
## DaysF111:CommunityD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Emmeans
Ssid.YII.emm<-emmeans(LM_Ssid.qpCR, ~Treatment*DaysF*Community)
#contrast(Ssid.YII.emm, "tukey")
# Pairwise comparisons
#pairs(Ssid.YII.emm) # same than contrast(Ssid.YII.emm, "tukey")
plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF), comparisons = TRUE) +
theme_bw() # Tukey comparission (do not trust CI to compare EMMs)
plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF|Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(Community~DaysF) +
theme_bw()
#CLD
SH_Multicomp<-cld(Ssid.YII.emm, by=NULL) # compact-letter display
SH_Multicomp<-SH_Multicomp[order(SH_Multicomp$DaysF,
SH_Multicomp$Community,
SH_Multicomp$Treatment), ]
#write.csv(SH_Multicomp, "Outputs/Ssid_SH_Multicomp.csv") # 2. Predict values:
pred_Ssid1 <- predict(LM_Ssid.qpCR, re.form = NA)
#3. Bootstrap CI:
Ss.boot1 <- bootMer(LM_Ssid.qpCR, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.boot1$t, 2, sd)
CI.lo_1 <- pred_Ssid1 - std.err*1.96
CI.hi_1 <- pred_Ssid1 + std.err*1.96
#Plot
Model_Ss_1b_plot<- ggplot(
qPCR.variables_3, aes(x = Days, y = logSH, colour = Treatment)) +
geom_line(aes(y = pred_Ssid1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
#scale_y_continuous(name=expression(~italic("Fv / Fm")),
# limits = c(0.5, 0.61),
# breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
#scale_x_continuous("Days in the experiment", limits = c(0, 78),
# breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1)+
# stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_1b_plot + facet_wrap(~Community)
- Days as continous
LM_Ssid.qpCR2<-lmer(logSH ~ Treatment * Days * Community + (1|Replicate) + (1|Genotype/Fragment),
data= qPCR.variables_2)
step(LM_Ssid.qpCR2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 -568.24 1180.5
## (1 | Replicate) 1 21 -568.24 1178.5 0.00 1 1
## (1 | Fragment:Genotype) 2 20 -568.24 1176.5 0.00 1 1
## (1 | Genotype) 0 19 -589.20 1216.4 41.93 1 9.461e-11
##
## <none>
## (1 | Replicate)
## (1 | Fragment:Genotype)
## (1 | Genotype) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:Days:Community 1 1.7016 0.4254 4 461.15 0.9054
## Days:Community 2 1.7500 0.8750 2 465.80 1.8642
## Treatment:Community 3 3.2309 0.8077 4 467.27 1.7133
## Community 0 22.9355 11.4678 2 364.82 24.1860
## Treatment:Days 0 3.0544 1.5272 2 471.54 3.2209
## Pr(>F)
## Treatment:Days:Community 0.4606
## Days:Community 0.1562
## Treatment:Community 0.1458
## Community 1.369e-10 ***
## Treatment:Days 0.0408 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days anova(LM_Ssid.qpCR2) ranova(LM_Ssid.qpCR2)LM_Ssid.qpCR2<-lmer(logSH ~Treatment + Days + Community + (1|Genotype) + Treatment:Days,
data= qPCR.variables_2)
step(LM_Ssid.qpCR2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 10 -542.81 1105.6
## (1 | Genotype) 0 9 -564.04 1146.1 42.471 1 7.173e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Community 0 22.9355 11.4678 2 364.82 24.1860 1.369e-10
## Treatment:Days 0 3.0544 1.5272 2 471.54 3.2209 0.0408
##
## Community ***
## Treatment:Days *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days anova(LM_Ssid.qpCR2) ranova(LM_Ssid.qpCR2) summary(LM_Ssid.qpCR2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
## Data: qPCR.variables_2
##
## REML criterion at convergence: 1085.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5286 -0.6117 -0.1087 0.6507 2.3774
##
## Random effects:
## Groups Name Variance Std.Dev.
## Genotype (Intercept) 0.6803 0.8248
## Residual 0.4741 0.6886
## Number of obs: 486, groups: Genotype, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.853179 0.360321 8.344461 -5.143 0.000773 ***
## TreatmentN -0.096926 0.136250 471.094970 -0.711 0.477201
## TreatmentN+P -0.266358 0.143608 471.050405 -1.855 0.064257 .
## Days -0.008168 0.001209 471.227691 -6.757 4.18e-11 ***
## CommunityC3 1.818818 0.283439 290.421477 6.417 5.64e-10 ***
## CommunityD 0.262047 0.169203 453.518949 1.549 0.122148
## TreatmentN:Days 0.003229 0.001701 471.924639 1.898 0.058301 .
## TreatmentN+P:Days 0.004293 0.001799 471.108284 2.386 0.017420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnN TrtN+P Days CmmnC3 CmmntD TrtN:D
## TreatmentN -0.183
## TreatmntN+P -0.185 0.480
## Days -0.248 0.591 0.564
## CommunityC3 -0.371 -0.026 0.005 0.053
## CommunityD -0.382 -0.009 0.011 0.047 0.581
## TrtmntN:Dys 0.149 -0.834 -0.400 -0.706 0.063 -0.018
## TrtmntN+P:D 0.141 -0.399 -0.839 -0.669 0.021 0.020 0.477 # 2. Predict values:
pred_Ssid1 <- predict(LM_Ssid.qpCR2,re.form = NA)
#3. Bootstrap CI:
Ss.boot1 <- bootMer(LM_Ssid.qpCR2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.boot1$t, 2, sd)
CI.lo_1 <- pred_Ssid1 - std.err*1.96
CI.hi_1 <- pred_Ssid1 + std.err*1.96
#Plot
Model_Ss_1b_plot<- ggplot(
qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
geom_line(aes(y = pred_Ssid1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
#scale_y_continuous(name=expression(~italic("Fv / Fm")),
# limits = c(0.5, 0.61),
# breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
#scale_x_continuous("Days in the experiment", limits = c(0, 78),
# breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1) +
# stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_1b_plot + facet_wrap(Genotype~Community)